;+
; NAME:
;       MAKE_Z_PLOT
;
; PURPOSE:
;       Creates the Z (vertical position) plots for the EMAF sample.
;
; CATEGORY:
;       EMAF plotting routine
;
; CALLING SEQUENCE:
;       MAKE_Z_PLOT
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       NONE
;
; KEYWORD PARAMETERS:
;       NONE
;
; OUTPUTS:
;       NONE  [Creates an .eps file]
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       NONE
;
; MODIFICATION HISTORY:
;
;       Created:  05/08/12, TPEB -- Initial Version.
;       Modified: 11/13/12, TPEB -- Moved to emaf/ directory and made
;                                   compliant.
;       Modified: 03/13/13, TPEB -- Updated to OMNI* framework.
;       Modified: 03/18/13, TPEB -- Added B&W version for published
;                                   paper to avoid paying $$$.
;
;-

PRO MAKE_Z_PLOT
  
  restore,'./local/BGPS_V1_emaf.sav',/ver
  restore,'./local/BGPS_V1_pvec.sav',/ver
  eind = where(constrain.emaf)
  pvec = pvec[eind]
  constrain = constrain[eind]

  ;; Restrict consideration to "well-constrained" sources...
  cind = where(constrain.post)
  emaf = emaf[cind]
  pvec = pvec[cind]
  constrain = constrain[cind]
  n = n_elements(emaf)
  
  
  ;; emaf.grs_match = abs(pvec.stat.duse[0] - emaf.grs_dist) LE 1.
  
  
  omni_lbd2rz,emaf.l,emaf.b,pvec.stat.duse[0],r,z
  pvec.stat.duse /= 1.d3
  
  
  ;; r = ruse
  ;; z = zuse
  
  
  ;;===========================================
  ;; Make the PS plot
  
  ;; myps,'./emaf_paper/plots/z_position_subs.eps'
  
  ;; ;; DIVIDE the group into 3 sections by Rgal
  rs = r[sort(r)]/1.d3
  ;; zs = z[sort(r)]
  
  cut = [0,long((findgen(3)+1.)/3.*n)-1]
  ;; print,rs[cut],format="(F0.2)"
  ;; print,cut,n_elements(zs)
  
  pm = ' '+cgSymbol('+-')+' '
  ;; xrange = [-150,150]
  
  ;; ;; hist = HISTOGRAM(z, binsize=5, loc=zout)
  ;; plothist,zs[cut[0]:cut[1]],zout,hist,psym=10,ytit='N per 5 pc bin',thick=3,$
  ;;          xtit='Z (Vertical Position) [pc]',xr=xrange,/xst,bin=5,charsize=1.0
  ;; yfit = mpfitpeak(zout, hist, A, NTERMS=3, ERROR=sqrt(float(hist)),PERROR=perr)
  ;; oplot,zout,yfit,color=cgColor('Black'),thick=3
  ;; txtstr = ['Fit centroid = '+string(A[1],pm,perr[1],format=$
  ;;                                    "(F+0.1,A0,F0.1)")+' pc',$
  ;;           'FWHM = '+string(A[2]*sqrt(8. * alog(2.)),format="(F0.1)")+' pc',$
  ;;           'Gauss Sig = '+string(A[2],pm,perr[2],format="(F0.1,A0,F0.1)")+' pc',$
  ;;           'R!dgal!n = '+string(rs[cut[0]],format="(F0.2)")+' - '+$
  ;;           string(rs[cut[1]],format="(F0.2,' kpc')")]
  
  ;; al_legend,txtstr,box=0,/top,/left,charsize=1.0,spacing=1.5
  
  
  ;; plothist,zs[cut[1]+1:cut[2]],zout,hist,psym=10,bin=5,charsize=1.0,/overplot,$
  ;;          color=cgColor('DeepPink'),thick=3
  ;; yfit = mpfitpeak(zout, hist, A, NTERMS=3, ERROR=sqrt(float(hist)),PERROR=perr)
  ;; oplot,zout,yfit,color=cgColor('DeepPink'),thick=3
  ;; txtstr = ['Fit centroid = '+string(A[1],pm,perr[1],format=$
  ;;                                    "(F+0.1,A0,F0.1)")+' pc',$
  ;;           'FWHM = '+string(A[2]*sqrt(8. * alog(2.)),format="(F0.1)")+' pc',$
  ;;           'Gauss Sig = '+string(A[2],pm,perr[2],format="(F0.1,A0,F0.1)")+' pc',$
  ;;           'R!dgal!n = '+string(rs[cut[1]+1],format="(F0.2)")+' - '+$
  ;;           string(rs[cut[2]],format="(F0.2,' kpc')")]
  
  ;; al_legend,txtstr,box=0,/top,/right,charsize=1.0,spacing=1.5,textcolor='DeepPink'
  
  
  ;; plothist,zs[cut[2]+1:cut[3]],zout,hist,psym=10,bin=5,charsize=1.0,/overplot,$
  ;;          color=cgColor('DodgerBlue'),thick=3
  ;; yfit = mpfitpeak(zout, hist, A, NTERMS=3, ERROR=sqrt(float(hist)),PERROR=perr)
  ;; oplot,zout,yfit,color=cgColor('DodgerBlue'),thick=3
  ;; txtstr = ['Fit centroid = '+string(A[1],pm,perr[1],format=$
  ;;                                    "(F+0.1,A0,F0.1)")+' pc',$
  ;;           'FWHM = '+string(A[2]*sqrt(8. * alog(2.)),format="(F0.1)")+' pc',$
  ;;           'Gauss Sig = '+string(A[2],pm,perr[2],format="(F0.1,A0,F0.1)")+' pc',$
  ;;           'R!dgal!n = '+string(rs[cut[2]+1],format="(F0.2)")+' - '+$
  ;;           string(rs[cut[3]],format="(F0.2,' kpc')")]
  
  ;; al_legend,txtstr,box=0,/bottom,/left,charsize=1.0,spacing=1.5,$
  ;;           textcolor='DodgerBlue'
  
  ;; axis,xaxis=0,xthick=3,xr=xrange,/xst,xtickformat='blank_axis'
  
  ;; myps,/done
  
  
  ;;===========================================
  ;; Make the PS plot
  
  myps,'./emaf_paper/plots/z_position.eps',xsize=3.1,ysize=3.1,/bw
  
  pm = ' '+cgSymbol('+-')+' '
  fwhm = sqrt(8.d*alog(2.d))
  ;; multiplot,[2,1],/doyaxis
  xrange=[-149.5,149.5]
  yrange=[0,80]
  cs = 0.7
  
  plothist,z,zout,hist,psym=10,xtit='z!dgal!n  [pc]',$
           ytit='N per 5 pc bin',thick=3,xr=xrange,/xst,bin=5,charsize=cs,$
           yr=yrange,/yst,ytickformat="(I0)",$
           /fill,fcolor='BLK5',color='BLK5',xmargin=[7,1],ymargin=[4,1]

  yfit = mpfitpeak(zout, hist, A, NTERMS=3, ERROR=sqrt(float(hist)),PERROR=perr)
  print,'HIST: ',m4_stat(hist)
  zc = findgen(1001)/5. - 100
  yc = gauss1(zc,[A[1],A[2],A[0]],/peak)
  cgOplot,zc,yc,thick=3
  ;; plothist,z,zout,hist,psym=10,xtit='Z (Vertical Position) [pc]',$
  ;;          ytit='N per 5 pc bin',thick=3,xr=xrange,/xst,bin=5,charsize=1.0,$
  ;;          title='Full Sample',yr=yrange,/yst,ytickformat="(I0)",/overplot
  txtstr = ['Fit centroid = '+string(A[1],pm,perr[1],format=$
                                     "(F+0.1,A0,F0.1)")+' pc',$
            'FWHM = '+string(A[2]*fwhm,pm,perr[2]*fwhm,format="(F0.1,A0,F0.1)")+' pc'];; ,$
            ;; 'Gauss Sig = '+string(A[2],pm,perr[2],format="(F0.1,A0,F0.1)")+' pc',$
            ;; 'R!dgal!n = '+string(rs[cut[0]],format="(F0.2)")+' - '+$
            ;; string(rs[cut[3]],format="(F0.2,' kpc')")]
  
  ;; al_legend,txtstr,box=0,/top,/left,charsize=cs*0.8,spacing=1.3
  FOR iii=0,n_elements(txtstr)-1 DO BEGIN
     cgText,alignment=0,orientation=0,color='black',charsize=cs*0.8,$
            -130,73.5-(iii*4),txtstr[iii]
     
  ENDFOR

  al_legend,['N = '+string(n_elements(emaf),format="(I0)")],box=0,/top,/right,$
            charsize=cs
  
  axis,xaxis=0,xthick=3,xr=xrange,/xst,xtickformat='blank_axis'
  
  myps,/done

  ;; multiplot
  
  
  ;; ;;===========================================
  ;; ;; Make the panel for GRS-YES sources only!
  
  
  ;; ind = WHERE(emaf.grs AND emaf.grs_match, nind)
  ;; print,nind
  ;; print,fix(total(emaf.grs))-nind
  
  
  
  ;; plothist,z[ind],zout,hist,psym=10,xtit='Z (Vertical Position) [pc]',$
  ;;          ;;ytit='N per 5 pc bin',$
  ;;          thick=3,xr=xrange,/xst,bin=5,charsize=1.0,$
  ;;          title='GRS Distance Matched',yr=yrange,/yst,$
  ;;          /fill,fcolor='BLK5',color='BLK5'


  ;; yfit = mpfitpeak(zout, hist, A, NTERMS=3, ERROR=sqrt(float(hist)),PERROR=perr)
  ;; yc = gauss1(zc,[A[1],A[2],A[0]],/peak)
  ;; cgOplot,zc,yc,thick=3
  ;; txtstr = ['Fit centroid = '+string(A[1],pm,perr[1],format=$
  ;;                                    "(F+0.1,A0,F0.1)")+' pc',$
  ;;           'FWHM = '+string(A[2]*fwhm,pm,perr[2]*fwhm,format="(F0.1,A0,F0.1)")+' pc',$
  ;;           'Gauss Sig = '+string(A[2],pm,perr[2],format="(F0.1,A0,F0.1)")+' pc',$
  ;;           'R!dgal!n = '+string(min(r[ind]/1.d3),format="(F0.2)")+' - '+$
  ;;           string(max(r[ind]/1.d3),format="(F0.2,' kpc')")]
  ;; ;; plothist,z[ind],zout,hist,psym=10,xtit='Z (Vertical Position) [pc]',$
  ;; ;;          ;;ytit='N per 5 pc bin',$
  ;; ;;          thick=3,xr=xrange,/xst,bin=5,charsize=1.0,$
  ;; ;;          title='GRS Distance Matched',yr=yrange,/yst,/overplot
  
  ;; al_legend,txtstr,box=0,/top,/left,charsize=0.8,spacing=1.3
  ;; al_legend,['N = '+string(nind,format="(I0)")],box=0,/top,/right,charsize=0.8
  
  ;; axis,xaxis=0,xthick=3,xr=xrange,/xst,xtickformat='blank_axis'
  
  ;; myps,/done,/mp
  
  
  
  ;;====================================================
  ;; Make plot for Z vs. \dsun
  myps,'./emaf_paper/plots/z_vs_dsun_bw.eps',xsize=3.1,ysize=3.1,/bw;ct=34,/reverse
  
  cs = 0.7
  
  plot,charsize=cs,pvec.stat.duse[0],z,psym=16,symsize=0.4,/nodata,$
         yr=[0,15],yst=4,$;; xmargin=[10,10],$
         xtit='d'+cgSymbol('sun')+'  [kpc]',/isotropic,ymargin=[4,0],$
       xmargin=[8,2]

  cgAxis,yaxis=0,ytit='z!dgal!n  [pc]',yr=[-150,100],/save,charsize=cs
  cgAxis,yaxis=1,/yst,ytickformat='blank_axis'
  
  dp = findgen(101)/100.*!x.crange[1]
  omni_lbd2rz,30,+0.5,dp*1.d3,rp1,zp1
  omni_lbd2rz,30,-0.5,dp*1.d3,rp2,zp2
  cgOplot,dp,zp1,color='blk5',thick=3,linestyle=2 ; Cyan in color version
  cgOplot,dp,zp2,color='blk5',thick=3,linestyle=2 ; Cyan in color version
  
  
  vline,/h,+60,color='blk5',linestyle=4,thick=3
  vline,/h,-60,color='blk5',linestyle=4,thick=3
  
  
  ;; col = fix(round( (emaf.l-7.5) / 57.5 * 255.) )
  ;; FOR jj=0,n_elements(col)-1 DO $
  ;;    cgOplot,duse[jj,0],zuse[jj],psym=16,symsize=0.4,color=col[jj]
  cgOplot,pvec.stat.duse[0],z,psym=16,symsize=0.35
  
  
  ;; cgColorbar,/right,/vert,title='GLON [deg]',range=[7.5,65],divi=0,$
  ;;            pos=[0.88,!y.window[0],0.92,!y.window[1]],charsize=cs
  
  myps,/done

END
